
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# DISCLAIMER AND GENERAL INFORMATION
#
# File: App_D_heterogeneous_treatment_effects.R
# Purpose: Analysis of heterogeneous treatment effects of main results
# Date: 10 July 2024
# Data: Survey data (pulled through 0_cbam_prep.R)
#
# Technical disclaimer:
# All analyses in R version 4.4.1 (2024-06-14 ucrt) -- "Race for Your Life"
# R Studio 2024.04.2 Build 764 ("Chocolate Cosmos" Release (e4392fc9, 2024-06-05) for Windows)
# Windows 10 Enterprise, 64-bit
# 12th Gen Intel(R) Core(TM) i7-1255U 1.70 GHz with 16GB RAM
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++



# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (A) Load data and packages ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Run source file to load data
source("0_cbam_prep.R")

# Load packages
library("tidyverse")
library("janitor")
library("ggpubr")
library("modelsummary")
options(modelsummary_format_numeric_latex = "plain")
library("marginaleffects")
library("xtable")
library("tinytable")


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (B) Heterogeneous treatment effects: Income ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

dta$income_clean <- ifelse(dta$income<14,dta$income,NA) 

dta %>%
  group_by(country) %>%
  summarise(median=median(income_clean, na.rm=TRUE))

dta$income_high <- ifelse(dta$income_clean>9 & dta$country=="germany", 1,
                   ifelse(dta$income_clean>8 & dta$country=="hungary", 1,
                   ifelse(dta$income_clean>6 & dta$country=="switzerland", 1,
                   ifelse(dta$income_clean>7 & dta$country=="uk", 1, 0))))


# Regression models: all conditions
ols_de0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$income_high==0,]) 
ols_de1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$income_high==1,]) 
ols_hu0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$income_high==0,]) 
ols_hu1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$income_high==1,]) 
ols_ch0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$income_high==0,]) 
ols_ch1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$income_high==1,]) 
ols_uk0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$income_high==0,]) 
ols_uk1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$income_high==1,]) 


order_full <- c("factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)4"="Treatment: Both")  


models.full <- list(
  "Germany: low income" = ols_de0,
  "Germany: high income" = ols_de1,
  "Hungary: low income" = ols_hu0,
  "Hungary: high income" = ols_hu1,
  "Switzerland: low income" = ols_ch0,
  "Switzerland: high income" = ols_ch1,
  "UK: low income" = ols_uk0,
  "UK: high income" = ols_uk1)

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- modelplot(rev(models.full), coef_map=order_full, coef_omit="Interc", conf_level=0.834, facet=TRUE) +
  labs(x="Coefficients",
       title="Treatment Effects across Countries",
       subtitle="DV: CBAM support (1-5 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic()
p1

# Export plot to PDF (Figure D.1)
ggexport(p1, filename="./hetero_income.pdf", width=5, height=7)




# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (C) Heterogeneous treatment effects: Education ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

dta %>%
  group_by(country) %>%
  summarise(median=median(edu_cat, na.rm=TRUE))

dta$edu_high <- ifelse(dta$edu_cat>2,1,0)

# Regression models: all conditions
ols_de0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$edu_high==0,]) 
ols_de1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$edu_high==1,]) 
ols_hu0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$edu_high==0,]) 
ols_hu1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$edu_high==1,]) 
ols_ch0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$edu_high==0,]) 
ols_ch1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$edu_high==1,]) 
ols_uk0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$edu_high==0,]) 
ols_uk1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$edu_high==1,]) 


order_full <- c("factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)4"="Treatment: Both")  


models.full <- list(
  "Germany: below median education" = ols_de0,
  "Germany: above median education" = ols_de1,
  "Hungary: below median education" = ols_hu0,
  "Hungary: above median education" = ols_hu1,
  "Switzerland: below median education" = ols_ch0,
  "Switzerland: above median education" = ols_ch1,
  "UK: below median education" = ols_uk0,
  "UK: above median education" = ols_uk1)

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- modelplot(rev(models.full), coef_map=order_full, coef_omit="Interc", conf_level=0.834, facet=TRUE) +
  labs(x="Coefficients",
       title="Treatment Effects across Countries",
       subtitle="DV: CBAM support (1-5 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic()
p1

# Export plot to PDF (Figure D.2)
ggexport(p1, filename="./hetero_education.pdf", width=5, height=7)



# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (D) Heterogeneous treatment effects: Public/private sector ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++


# Regression models: all conditions
ols_de0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$sector==1,]) 
ols_de1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & (dta$sector==2 | dta$sector==3),]) 
ols_hu0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$sector==1,]) 
ols_hu1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & (dta$sector==2 | dta$sector==3),]) 
ols_ch0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$sector==1,]) 
ols_ch1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & (dta$sector==2 | dta$sector==3),]) 
ols_uk0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$sector==1,]) 
ols_uk1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & (dta$sector==2 | dta$sector==3),]) 


order_full <- c("factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)4"="Treatment: Both")  


models.full <- list(
  "Germany: public sector" = ols_de0,
  "Germany: private sector" = ols_de1,
  "Hungary: public sector" = ols_hu0,
  "Hungary: private sector" = ols_hu1,
  "Switzerland: public sector" = ols_ch0,
  "Switzerland: private sector" = ols_ch1,
  "UK: public sector" = ols_uk0,
  "UK: private sector" = ols_uk1)

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- modelplot(rev(models.full), coef_map=order_full, coef_omit="Interc", conf_level=0.834, facet=TRUE) +
  labs(x="Coefficients",
       title="Treatment Effects across Countries",
       subtitle="DV: CBAM support (1-5 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic()
p1

# Export plot to PDF (Figure D.3)
ggexport(p1, filename="./hetero_sector.pdf", width=5, height=7)



# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (E) Heterogeneous treatment effects: Retired ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Regression models: all conditions
ols_de0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$retired==0,]) 
ols_de1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$retired==1,]) 
ols_hu0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$retired==0,]) 
ols_hu1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$retired==1,]) 
ols_ch0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$retired==0,]) 
ols_ch1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$retired==1,]) 
ols_uk0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$retired==0,]) 
ols_uk1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$retired==1,]) 


order_full <- c("factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)4"="Treatment: Both")  


models.full <- list(
  "Germany: working age" = ols_de0,
  "Germany: retired" = ols_de1,
  "Hungary: working age" = ols_hu0,
  "Hungary: retired" = ols_hu1,
  "Switzerland: working age" = ols_ch0,
  "Switzerland: retired" = ols_ch1,
  "UK: working age" = ols_uk0,
  "UK: retired" = ols_uk1)

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- modelplot(rev(models.full), coef_map=order_full, coef_omit="Interc", conf_level=0.834, facet=TRUE) +
  labs(x="Coefficients",
       title="Treatment Effects across Countries",
       subtitle="DV: CBAM support (1-5 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic()
p1

# Export plot to PDF (Figure D.4)
ggexport(p1, filename="./hetero_retired.pdf", width=5, height=7)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (F) Heterogeneous treatment effects: Unemployed ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Regression models: all conditions
ols_de0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$wowork==0,]) 
ols_de1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$wowork==1,]) 
ols_hu0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$wowork==0,]) 
ols_hu1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$wowork==1,]) 
ols_ch0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$wowork==0,]) 
ols_ch1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$wowork==1,]) 
ols_uk0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$wowork==0,]) 
ols_uk1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$wowork==1,]) 


order_full <- c("factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)4"="Treatment: Both")  


models.full <- list(
  "Germany: working" = ols_de0,
  "Germany: umployed" = ols_de1,
  "Hungary: working" = ols_hu0,
  "Hungary: unemployed" = ols_hu1,
  "Switzerland: working" = ols_ch0,
  "Switzerland: unemployed" = ols_ch1,
  "UK: working" = ols_uk0,
  "UK: unemployed" = ols_uk1)

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- modelplot(rev(models.full), coef_map=order_full, coef_omit="Interc", conf_level=0.834, facet=TRUE) +
  labs(x="Coefficients",
       title="Treatment Effects across Countries",
       subtitle="DV: CBAM support (1-5 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic()
p1

# Export plot to PDF (Figure D.5)
ggexport(p1, filename="./hetero_unemployed.pdf", width=5, height=7)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (G) Heterogeneous treatment effects: Climate concern ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

dta %>%
  group_by(country) %>%
  summarise(median=median(cc_concern, na.rm=TRUE))

dta$cc_high <- ifelse(dta$cc_concern>3,1,0)


# Regression models: all conditions
ols_de0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$cc_high==0,]) 
ols_de1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$cc_high==1,]) 
ols_hu0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$cc_high==0,]) 
ols_hu1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$cc_high==1,]) 
ols_ch0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$cc_high==0,]) 
ols_ch1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$cc_high==1,]) 
ols_uk0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$cc_high==0,]) 
ols_uk1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$cc_high==1,]) 


order_full <- c("factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)4"="Treatment: Both")  

models.full <- list(
  "Germany: low concern" = ols_de0,
  "Germany: high concern" = ols_de1,
  "Hungary: low concern" = ols_hu0,
  "Hungary: high concern" = ols_hu1,
  "Switzerland: low concern" = ols_ch0,
  "Switzerland: high concern" = ols_ch1,
  "UK: low concern" = ols_uk0,
  "UK: high concern" = ols_uk1)

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- modelplot(rev(models.full), coef_map=order_full, coef_omit="Interc", conf_level=0.834, facet=TRUE) +
  labs(x="Coefficients",
       title="Treatment Effects across Countries",
       subtitle="DV: CBAM support (1-5 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic()
p1

# Export plot to PDF (Figure D.6)
ggexport(p1, filename="./hetero_ccconcern.pdf", width=5, height=7)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (H) Heterogeneous treatment effects: Left-right ideology ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Regression models: all conditions
ols_de0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$lr_cat==1,]) 
ols_de1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$lr_cat==3,]) 
ols_hu0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$lr_cat==1,]) 
ols_hu1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$lr_cat==3,]) 
ols_ch0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$lr_cat==1,]) 
ols_ch1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$lr_cat==3,]) 
ols_uk0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$lr_cat==1,]) 
ols_uk1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$lr_cat==3,]) 


order_full <- c("factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)4"="Treatment: Both")  


models.full <- list(
  "Germany: left" = ols_de0,
  "Germany: right" = ols_de1,
  "Hungary: left" = ols_hu0,
  "Hungary: right" = ols_hu1,
  "Switzerland: left" = ols_ch0,
  "Switzerland: right" = ols_ch1,
  "UK: left" = ols_uk0,
  "UK: right" = ols_uk1)

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- modelplot(rev(models.full), coef_map=order_full, coef_omit="Interc", conf_level=0.834, facet=TRUE) +
  labs(x="Coefficients",
       title="Treatment Effects across Countries",
       subtitle="DV: CBAM support (1-5 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic()
p1

# Export plot to PDF (Figure D.7)
ggexport(p1, filename="./hetero_leftright.pdf", width=5, height=7)

# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (I) Heterogeneous treatment effects: Urban-rural ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Regression models: all conditions
ols_de0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$population_density==3,]) 
ols_de1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$population_density<3,]) 
ols_hu0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$population_density==3,]) 
ols_hu1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$population_density<3,]) 
ols_ch0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$population_density==3,]) 
ols_ch1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$population_density<3,]) 
ols_uk0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$population_density==3,]) 
ols_uk1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$population_density<3,]) 


order_full <- c("factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)4"="Treatment: Both")  


models.full <- list(
  "Germany: rural" = ols_de0,
  "Germany: urban" = ols_de1,
  "Hungary: rural" = ols_hu0,
  "Hungary: urban" = ols_hu1,
  "Switzerland: rural" = ols_ch0,
  "Switzerland: urban" = ols_ch1,
  "UK: rural" = ols_uk0,
  "UK: urban" = ols_uk1)

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- modelplot(rev(models.full), coef_map=order_full, coef_omit="Interc", conf_level=0.834, facet=TRUE) +
  labs(x="Coefficients",
       title="Treatment Effects across Countries",
       subtitle="DV: CBAM support (1-5 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic()
p1

# Export plot to PDF (Figure D.8)
ggexport(p1, filename="./hetero_ruralurban.pdf", width=5, height=7)


# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# (J) Heterogeneous treatment effects: Regional differences ----
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

# Germany: East and West
dta$de_east <- ifelse(dta$country=="germany" & (dta$region=="Brandenburg" |
                                                dta$region=="Mecklenburg-Vorpommern" |
                                                dta$region=="Sachsen" |
                                                dta$region=="Sachsen-Anhalt" |
                                                dta$region=="Thüringen"),1,
                      ifelse(dta$country=="germany",0,NA))

# Hungary: Eastern, Western, and Central region
dta$hu_west <- ifelse(dta$country=="hungary" & dta$region=="Dunántúl", 1,
                      ifelse(dta$country=="hungary",0,NA))

dta$hu_centre <- ifelse(dta$country=="hungary" & dta$region=="Közép-Magyarország", 1,
                      ifelse(dta$country=="hungary",0,NA))

dta$hu_east <- ifelse(dta$country=="hungary" & dta$region=="Alföld és Észak", 1,
                        ifelse(dta$country=="hungary",0,NA))

# Switzerland: East and West
dta$ch_west <- ifelse(dta$country=="switzerland" & (dta$county=="Genève" |
                                                  dta$county=="(le) Jura" |
                                                  dta$county=="Fribourg" |
                                                  dta$county=="Fribourg / Freiburg" |
                                                  dta$county=="Neuchâtel" |
                                                  dta$county=="Valais / Wallis" |
                                                  dta$county=="Vaud"), 1,
                      ifelse(dta$country=="switzerland",0,NA))


# UK: Northern and Southern region
dta$uk_north <- ifelse(dta$country=="uk" & (dta$county=="Scotland" |
                                            dta$county=="North West" |
                                            dta$county=="North East" |
                                            dta$county=="Yorkshire and Humberside") , 1,
                      ifelse(dta$country=="uk",0,NA))


# Regression models: all conditions
ols_de0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$de_east==1,]) 
ols_de1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="germany" & dta$de_east==0,]) 
ols_hu1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$hu_east==1,]) 
ols_hu2 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$hu_centre==1,])
ols_hu3 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="hungary" & dta$hu_west==1,])
ols_ch0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$ch_west==0,]) 
ols_ch1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="switzerland" & dta$ch_west==1,]) 
ols_uk0 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$uk_north==1,]) 
ols_uk1 <- lm(cbam_support~factor(tr_main), data=dta[dta$country=="uk" & dta$uk_north==0,]) 

order_full <- c("factor(tr_main)2"="Treatment: Jobs",
                "factor(tr_main)3"="Treatment: Prices",
                "factor(tr_main)4"="Treatment: Both")  

models.full <- list(
  "Germany: East" = ols_de0,
  "Germany: West" = ols_de1,
  "Hungary: East" = ols_hu1,
  "Hungary: Central" = ols_hu2,  
  "Hungary: West" = ols_hu2,
  "Switzerland: East" = ols_ch0,
  "Switzerland: West" = ols_ch1,
  "UK: North" = ols_uk0,
  "UK: South" = ols_uk1)

# Plot coefficients

# Set color-blind friendly color palatte
# http://www.cookbook-r.com/Graphs/Colors_(ggplot2)/#a-colorblind-friendly-palette
cols <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")

p1 <- modelplot(rev(models.full), coef_map=order_full, coef_omit="Interc", conf_level=0.834, facet=TRUE) +
  labs(x="Coefficients",
       title="Treatment Effects across Countries",
       subtitle="DV: CBAM support (1-5 scale)") +
  geom_vline(xintercept=0, color='black', linetype = "dashed") +
  theme_classic()
p1

# Export plot to PDF (Figure D.9)
ggexport(p1, filename="./hetero_region.pdf", width=5, height=7)






# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#                         END OF FILE
# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
